#include "Matrix.hpp"

Matrix::Matrix(unsigned int rows, unsigned int cols) :
	_rows(rows), _cols(cols)
{
	createMatrix();
	// inicializa la matriz en cero
	__MATRIX_EACH(i,j) _matrix[i][j] = 0;
}

Matrix::Matrix(unsigned int rows, unsigned int cols, const double* arr) :
	_rows(rows), _cols(cols)
{
	createMatrix();
	__MATRIX_EACH(i,j) _matrix[i][j] = arr[i*cols+j];
}

Matrix::Matrix(unsigned int rows, unsigned int cols, const double** arr) :
	_rows(rows), _cols(cols)
{
	createMatrix();
	__MATRIX_EACH(i,j) _matrix[i][j] = arr[i][j];
}

Matrix::Matrix(const Matrix& m)
{
	_rows = m._rows;
	_cols = m._cols;
	
	_matrix = new double*[_rows];
	__MATRIX_EACHROW(i) _matrix[i] = new double[_cols];

	*this = m;
}

Matrix::~Matrix()
{
	// libera la memoria de la matriz
	__MATRIX_EACHROW(i) delete[] _matrix[i];
	delete[] _matrix;
}

void Matrix::createMatrix()
{
	if( _rows == 0 || _cols == 0 )
		throw std::domain_error( (char*)"El número de filas o columnas es cero." );

	// crea la matriz en memoria
	_matrix = new double*[_rows];
	__MATRIX_EACHROW(i) _matrix[i] = new double[_cols];
}

// cantidad de filas
unsigned int Matrix::rows() const
{
	return _rows;
}

// cantidad de columnas
unsigned int Matrix::cols() const
{
	return _cols;
}

Matrix Matrix::col(unsigned int j) const
{
	Matrix v(_rows,1);
	__MATRIX_EACHROW(i) v(i,1) = _matrix[i][j];
	return v;
}

Matrix Matrix::row(unsigned int i) const
{
	Matrix v(1,_cols);
	__MATRIX_EACHCOL(j) v(1,j) = _matrix[i][j];
	return v;
}

Matrix Matrix::minor(unsigned int i0, unsigned int j0) const
{
	Matrix result( _rows-1, _cols-1 );
	unsigned int poffset=0, qoffset=0;

	__MATRIX_EACHROW_M(result,p)
	{
		__MATRIX_EACHCOL_M(result,q)
		{
			if( p >= i0 )
				poffset = 1;
			if( q < j0 )
				qoffset = 0;
			else
				qoffset = 1;
			result._matrix[p][q] = _matrix[p+poffset][q+qoffset];
		}
	}

	return result;
}

// multiplicación de matrices
Matrix Matrix::operator* ( const Matrix& b ) const
{
	if( this->_cols != b._rows )
		throw std::logic_error( (char*)"Multiplicación de matrices incompatibles" );

	Matrix result( this->_rows, b._cols );

	__MATRIX_EACH_M(result,i,j)
	{
		double sum = 0;
		__MATRIX_EACHROW(k) sum += _matrix[i][k] * b._matrix[k][j];
		result._matrix[i][j] = sum;
	}

	return result;
}

// multiplicación por un escalar
Matrix Matrix::operator* ( const double k ) const
{
	Matrix result( this->_rows, this->_cols );

	__MATRIX_EACH_M(result,i,j) result._matrix[i][j] = this->_matrix[i][j]*k;

	return result;	
}

// suma
Matrix Matrix::operator+ ( const Matrix& b ) const
{
	if( this->_cols != b._cols || this->_rows != b._rows )
		throw std::logic_error( (char*)"Suma de matrices incompatibles" );

	Matrix result( this->_rows, this->_cols );

	__MATRIX_EACH_M(result,i,j) result._matrix[i][j] = this->_matrix[i][j]+b._matrix[i][j];

	return result;	
}

// resta
Matrix Matrix::operator- ( const Matrix& b ) const
{
	if( this->_cols != b._cols || this->_rows != b._rows )
		throw std::logic_error( (char*)"Resta de matrices incompatibles" );

	Matrix result( this->_rows, this->_cols );

	__MATRIX_EACH_M(result,i,j) result._matrix[i][j] = this->_matrix[i][j]-b._matrix[i][j];

	return result;	
}

// asignación
Matrix& Matrix::operator= ( const Matrix& b )
{
	if( this->_cols != b._cols || this->_rows != b._rows )
		throw std::logic_error( (char*)"Asignación de matrices incompatibles" );

	__MATRIX_EACH(i,j) this->_matrix[i][j] = b._matrix[i][j];

	return *this;
}

Matrix& Matrix::operator+= ( const Matrix& b )
{
	if( this->_cols != b._cols || this->_rows != b._rows )
		throw std::logic_error( (char*)"Suma de matrices incompatibles" );

	__MATRIX_EACH(i,j) this->_matrix[i][j] =  this->_matrix[i][j]+b._matrix[i][j];

	return *this;
}

Matrix& Matrix::operator-= ( const Matrix& b )
{
	if( this->_cols != b._cols || this->_rows != b._rows )
		throw std::logic_error( (char*)"Resta de matrices incompatibles" );

	__MATRIX_EACH(i,j) this->_matrix[i][j] =  this->_matrix[i][j]-b._matrix[i][j];

	return *this;
}

Matrix& Matrix::operator*= ( const double k )
{
	__MATRIX_EACH(i,j) this->_matrix[i][j] =  this->_matrix[i][j]*k;
	return *this;
}

Matrix& Matrix::operator/= ( const double k )
{
	__MATRIX_EACH(i,j) this->_matrix[i][j] = this->_matrix[i][j]/k;
	return *this;
}

Matrix& Matrix::operator*= ( const Matrix& b )
{
	if( this->_cols != b._rows )
		throw std::logic_error( (char*)"Multiplicación de matrices incompatibles" );

	__MATRIX_EACH(i,j)
	{
		double sum = 0;
		__MATRIX_EACHROW(k) sum += this->_matrix[i][k] * b._matrix[k][j];
		_matrix[i][j] = sum;
	}

	return *this;
}

double& Matrix::operator()( unsigned int i, unsigned int j )
{
	return _matrix[i-1][j-1];
}

const double Matrix::operator()( unsigned int i, unsigned int j ) const
{
	return _matrix[i-1][j-1];
}

// comparación
bool Matrix::operator== ( const Matrix& b ) const
{
	if( this->_cols != b._cols || this->_rows != b._rows )
		throw std::logic_error( (char*)"Comparación de matrices incompatibles" );

	bool eq = true;

	__MATRIX_EACH(i,j) eq &= (this->_matrix[i][j] == b._matrix[i][j]);

	return eq;
}

double Matrix::FrobeniusNorm() const
{
	double val = 0;
	__MATRIX_EACH(i,j) val += _matrix[i][j] * _matrix[i][j];
	return std::sqrt(val);
}

// representa a la matriz en un ostream
std::ostream& Matrix::print( std::ostream& out ) const
{
	__MATRIX_EACHROW(i)
	{
		__MATRIX_EACHCOL(j) out << _matrix[i][j] << ' ';
		out << std::endl;
	}
	return out;
}

Matrix Matrix::transposed() const
{
	Matrix result( _cols, _rows );

	__MATRIX_EACH(i,j) result._matrix[j][i] = _matrix[i][j];

	return result;
}

void Matrix::swapRows( unsigned int i1, unsigned int i2 )
{
	__MATRIX_EACHCOL(j)
	{
		double tmp = _matrix[i1][j];
		_matrix[i1][j] = _matrix[i2][j];
		_matrix[i2][j] = tmp;
	}
}

void Matrix::swapCols( unsigned int j1, unsigned int j2 )
{
	__MATRIX_EACHROW(i)
	{
		double tmp = _matrix[i][j1];
		_matrix[i][j1] = _matrix[i][j2];
		_matrix[i][j2] = tmp;
	}
}

void Matrix::substractRows( unsigned int i1, unsigned int i2, double k )
{
	__MATRIX_EACHCOL(j) _matrix[i1][j] -= k*_matrix[i2][j];
}

bool Matrix::isSquare() const
{
	return (_cols == _rows);
}

double Matrix::absmax( unsigned int& i0, unsigned int& j0 ) const
{
	unsigned int ii0 = i0, ij0 = j0;
	double max = std::abs(_matrix[i0][j0]);
	for( unsigned int i = ii0 ; i < _rows ; i++ )
		for( unsigned int j = ij0 ; j < _cols ; j++ )
		{
			double tmp = std::abs(_matrix[i][j]);
			if( tmp > max )
			{
				i0  = i;
				j0  = j;
				max = tmp;
			}
		}
	return max;
}

Matrix operator* ( const double k, const Matrix& m )
{
	return m*k;
}

std::istream& Matrix::read( std::istream& in )
{
	__MATRIX_EACH(i,j) in >> _matrix[i][j];
	return in;
}

std::ostream& operator<< ( std::ostream& out, const Matrix& m )
{
	return m.print(out);
}

std::istream& operator>> ( std::istream& in, Matrix& m )
{
	return m.read( in );
}

